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Abstract 

A  new  approach  to  frequency  domain  design  of  robust  controllers 
for  distributed  parameter  systems  is  presented.  The  central  idea  is 
to  use  techniques  that  were  developed  for  the  solution  of  the  Corona 
Problem,  for  the  solution  of  both  the  Bezout  equation  and  an  auxil¬ 
iary  equation  that  arises  form  the  Nehari  interpolation  problem.  An 
algebraic  reformulation  of  these  equations  allows  the  solution  to  be 
computed  from  the  solution  of  an  inhomogeneous  Cauchy  Riemann 
equation  with  a  Carleson  measure  as  the  inhomogeneous  term.  The 
theory  is  applied  to  a  single  input  single  output  system  with  delay  to 
yield  the  transfer  function  of  a  stabilizing  controller  with  guaranteed 
H°°  stability  margin.  Finally  the  framework  is  extended  to  handle 
multi-input  multi-output  systems. 
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1  Introduction 


This  paper  presents  a  new  approach  to  H00  controller  design  for  a  general 
class  of  linear  systems  that  are  described  by  irrational  transfer  functions. 
The  approach  is  based  on  a  computational  method  for  solving  linear  Dio- 
phantine  equations  in  algebras  over  the  ring  of  H00  functions.  The  term 
Diophantine  here  is  used  in  analogy  with  Diophantine  equations  over  the 
ring  of  integers  to  refer  to  algebraic  equations  over  more  general  rings,  in 
this  case  rings  constructed  from  function  spaces.  The  term  linear  is  used  to 
indicate  that  the  equations  that  are  being  considered  do  not  involve  products 
of  the  unknown  variable,  in  other  words  they  have  the  same  form  as  linear 
algebraic  equations  over  fields.  Two  such  equations  occur  in  the  design  of 
H00  controllers  for  linear  systems:  the  first  is  the  Bezout  identity  that  leads 
to  the  parameterization  of  all  stabilizing  controllers,  and  the  second  is  an 
equation  that  arises  from  the  Nehari  problem.  Bounded  solutions  to  the 
second  equation  provide  controllers  with  guaranteed  robustness  in  the  sense 
of  H°°  control  theory. 

The  extension  of  H°°  design  techniques  to  distributed  parameter  systems 
has  been  an  area  of  active  research  since  the  1980’s.  In  the  report  [Cur92] 
Curtain  compares  five  approaches  that  were  developed  during  this  period. 
The  approach  that  has  enjoyed  the  greatest  practical  success  is  based  on  the 
work  of  Glover,  Curtain  and  Partington  [GCP88],  in  this  paper  the  authors 
show  that  if  a  distributed  system  has  a  Hankel  operator  of  nuclear  type, 
then  the  system  may  be  approximated  by  a  finite  dimensional  system  with 
a  guaranteed  bound  on  the  L00  norm  of  the  error.  This  result  means  that 
a  finite  dimensional  controller  designed  for  the  approximate  plant  will  also 
control  the  infinite  dimensional  plant  with  guaranteed  (though  reduced)  H00 
stability  margins.  For  linear  plants  that  are  associated  with  a  nuclear  Hankel 
operator,  a  class  that  includes  many  cases  of  practical  interest,  this  method 
provides  a  practical  solution  to  the  problem  of  controller  design. 

In  spite  of  this  result,  research  has  continued  on  ways  to  attack  the 
problem  of  controller  design  directly,  without  the  initial  finite  dimensional 
approximation  to  the  plant.  As  in  the  case  of  the  theory  that  was  developed 
for  rational  plants,  the  work  on  infinite  dimensional  systems  can  be  divided 
between  approaches  that  consider  the  plant  as  a  linear  state-space  system 
with  an  infinite  dimensional  state  space,  and  approaches  that  consider  the 
plant  as  a  transfer  function  that  acts  multiplicatively  in  the  frequency  do¬ 
main  on  the  transform  of  the  signal  space.  An  example  of  the  first  approach 
is  found  in  the  work  of  van  Kuelen  [vK93]  in  which  the  author  presents  a 
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pair  of  infinite  dimensional  Riccati  equations  with  solutions  that  are  the  ba¬ 
sis  for  a  state  space  representation  of  a  controller.  Examples  of  the  second 
approach  are  found  in  the  work  of  Foias,  Tannenbaum,  Ozbay  and  Smith 
[FT88a],  [FT88bj,  [OST93],  [093],  and  in  the  work  of  Dym,  Geogiou  and 
Smith  [GS92],  [DGS93].  The  two  sets  of  authors  design  H°°  controllers  using 
methods  based  on  the  theory  of  skew  Toeplitz  operators  that  is  presented  by 
Bercovici,  Foias  and  Tannenbaum  in  [BFT88].  The  authors  Flamm,  Mitter 
and  Yang  independently  follow  a  similar  approach  in  [FM87]  and  [FY94], 

The  method  of  controller  design  presented  in  this  paper  takes  the  fre¬ 
quency  domain  approach.  It  extends  the  method  presented  by  Francis  in 
[Fra87]  for  systems  with  rational  transfer  functions  to  systems  with  irra¬ 
tional  transfer  functions.  The  method  starts  with  a  co-prime  factorization 
of  the  plant,  solves  the  associated  Bezout  Equation,  and  from  this  computes 
a  Youla  parameterization  of  all  stabilizing  controllers.  An  associated  Ne- 
hari  problem  is  solved  approximately  and  the  solution  provides  a  controller 
which,  although  not  optimal,  has  guaranteed  robustness  in  an  H°°  sense. 
The  result  of  the  controller  computation  is  a  numerical  approximation  to 
the  controller’s  transfer  functions  rather  than  a  closed  form  expression  com¬ 
posed  of  elementary  functions. 

The  main  contribution  of  this  paper  is  the  presentation  of  a  new  method 
for  solving  both  the  Bezout  equation  that  leads  to  the  Youla  parameteriza¬ 
tion  and  the  equation  that  is  derived  from  the  Nehari  Problem.  This  method 
is  based  on  a  constructive  proof  of  the  corona  theorem  presented  by  Gar¬ 
nett  in  [Gar81].  The  constructive  part  involves  an  explicit  solution  to  the 
inhomogeneous  Cauchy  Riemann  equation  db/dz  —  /r,  where  the  inhomoge¬ 
neous  part  /r  is  a  Carleson  measure,  and  the  solution  b  is  a  distribution  with 
support  in  the  right  half  plane  and  an  L°°  boundary  value.  The  terminology 
is  explained  in  detail  in  Section  3.  The  method  of  solution  of  this  first  order 
partial  differential  equation,  which  is  due  to  Jones  [Jon80]  [Jon83],  combines 
a  careful  decomposition  of  the  measure  /r  with  a  Green’s  function  method. 

It  seems  most  logical  to  compare  the  results  achieved  here  with  those 
achieved  by  the  methods  developed  in  [OST93],  and  to  this  end  the  example 
presented  by  Enns,  Ozbay  and  Tannenbaum  in  [EOT92]  is  reworked  in  this 
paper  with  the  techniques  developed  here.  The  steps  in  the  computation  of 
a  compensator  transfer  function  are  presented  in  detail,  and  an  open-loop 
Nyquist  plot  for  the  compensated  system  is  presented  along  with  magnitude 
plots  of  appropriate  sensitivity  transfer  functions  for  the  closed  loop  system. 
Enns  et  al.  are  able  to  compute  an  optimal  solution  for  this  problem,  so  a 
comparison  with  the  results  in  [EOT92]  provides  a  good  indication  of  the 
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distance  between  the  solution  computed  here  and  an  optimal  solution. 

The  paper  is  divided  into  six  sections  of  which  this  introduction  is  the 
first.  In  the  second  section  the  example  from  [EOT92]  is  presented  and  an 
outline  of  the  method  to  be  employed  for  its  solution  is  given,  in  particular 
this  section  shows  how  linear  Diophantine  equations  occur  in  the  H00  design 
problem.  The  third  section  provides  the  theory  that  is  used  to  solve  the 
linear  Diophantine  equations.  This  is  necessarily  the  most  technical  sec¬ 
tion  of  the  paper  and  includes  a  number  of  ideas  that  are  not  commonly 
met  in  the  control  literature.  The  fourth  section  shows  how  the  algorithms 
that  are  described  in  section  three  are  applied  to  the  example  introduced 
in  section  two.  Included  in  this  section  is  a  description  of  the  structure  of 
the  software  that  was  developed  to  implement  the  algorithms  and  the  re¬ 
sults  of  the  controller  design.  Section  five  contains  the  theory  needed  to 
extend  the  methods  presented  in  the  previous  sections  to  general  multiple- 
input  multiple-output  systems,  and  the  final  section  contains  provides  some 
conclusions  and  a  discussion  of  future  directions 
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2  An  Example 

A  concrete  example  is  introduced  in  this  section  in  order  to  establish  a  design 
method  and  provide  motivation  for  the  techniques  that  are  introduced  in  the 
following  section.  The  example  uses  a  simplified  model  of  the  pitch-axis  fast 
dynamics  of  an  unstable  aircraft  taken  from  a  paper  by  Enns,  Ozbay  and 
Tannenbaum  [EOT92].  The  plant,  which  is  given  by  the  transfer  function1 


F(z) 


e 


—  TZ 


(TZ  —  1  ’ 


(1) 


has  an  unstable  pole  at  z  =  1  /a,  and  a  delay  of  r  seconds,  as  such  it  is 
amongst  the  simplest  unstable  infinite  dimensional  systems.  The  design 
objective  is  to  produce  a  linear  feedback  controller  with  the  configuration 
illustrated  in  Figure  1  that  both  robustly  stabilizes  the  plant  and  maintains 
low  low-frequency  sensitivity. 


Figure  1:  Feedback  Controller 


The  design  proceeds  by  the  standard  procedure  [Fra87]  of  reformulating 
the  problem  as  the  minimization  of  the  L°°  norm  of  an  affine  expression  over 
an  H00  parameter.  Let 


Fi(z)  =  e~rz  /  (crz  +  1)  (2) 

Fi{z)  =  (crz  —  1)/ (crz  +  1) 

lThe  symbol  z  is  used  throughout  the  paper  to  denote  an  independent  complex  variable 
in,  for  example,  transfer  functions.  While  this  notation  is  at  odds  with  the  usual  choice 
of  s  for  the  independent  variable  in  the  Laplace  transform,  it  is  standard  in  complex 
analysis,  and  it’s  adoption  here  offers  the  benefit  of  maintaining  consistency  between  the 
notation  in  Section  3  and  the  notation  in  the  literature,  as  well  as  an  improvement  in  the 
self-consistency  of  the  notation. 
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then  F(z)  =  F\{z)/F2{z)  and  infRez>0(|i7i(,2)|  +  |F2(z)|)  >  0  so  F\  and  F2 
form  a  co-prime  factorization  of  F.  Let  X\  and  X2  be  functions  in  H°°  that 
solve  the  Bezout  equation2 

F1(z)X1(z)+F2(z)X2(z)  =  1  (3) 

then  all  stabilizing  controllers  may  be  expressed  in  terms  of  an  H00  parameter 
Q  by  the  bilinear  function 

X.jz)  +  F2(z)Q(z) 

X2(z)  -  Fi{z)Q{z) ' 

This  parameterization  was  first  given  by  Youla  et  al.  [YBJ76]  for  the  case  of 
rational  transfer  functions  and  was  extended  to  irrational  transfer  functions 
by  Baras  in  [Bar80].  Smith  gives  a  proof  of  the  existence  of  strongly  co-prime 
factorizations  for  stabilizable  plants  in  [Smi89]. 

A  key  step  in  the  parameterization  is  the  solution  of  the  Bezout  equation. 
If  the  plant  has  rational  co-prime  factors,  then  the  equation  may  be  solved 
by  algorithms  that  exploit  the  Euclidean  domain  structure  of  the  ring  of 
polynomials  such  as  that  given  by  Kailath  in  Chapter  2  of  [Kai80],  or  by 
algorithms  that  rely  on  state  space  techniques  such  as  that  used  by  Francis  in 
[Fra87].  Neither  technique  will  handle  irrational  transfer  functions  however, 
and  people  working  with  irrational  transfer  functions  have  been  restricted  to 
systems  that  are  sufficiently  simple  that  the  solution  to  the  Bezout  equation 
may  be  found  by  inspection.  This  is  the  case  in  the  example  presented  in 
Equation  (1);  the  authors  of  [EOT92]  give  the  solution 

eT/cr 

(az  +  1)  -  2 eT^e~TZ 
(az  -  1) 

The  lack  of  a  good  method  for  computing  solutions  to  the  Bezout  equation 
has  hindered  the  application  of  frequency  domain  methods  to  controller 
design  for  plants  with  irrational  transfer  functions. 

Returning  to  the  example,  and  in  line  with  the  H00  design  methodology, 
the  design  goals  of  robust  stabilization  and  low  low-frequency  sensitivity  are 
reformulated  as  the  requirement  that  the  controller  should  be  chosen  so  that 
the  weighted  norms  of  a  pair  of  transfer  functions  should  be  minimized.  The 
sensitivity  is  given  by  the  closed  loop  transfer  function  (1  +  FG)-1  and  the 

2  The  symbol  1  is  used  to  represent  the  constant  function  l(z)  =  1. 


Xi(z)  = 
X2(z)  = 
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robust  stabilization  requirement  is  interpreted  as  meaning  that  the  controller 
should  stabilize  all  plants  with  transfer  functions  within  an  L°°  neighbor¬ 
hood  of  the  nominal  plant.  An  argument  based  on  the  Nyquist  theorem 
(Theorem  1  of  [CD82])  assures  a  stable  neighborhood  provided  the  trans¬ 
fer  function  G(1  +  FG)~ 1  is  bounded  above.  Both  objectives  are  satisfied 
with  the  choice  of  a  feedback  controller  G  that  minimizes  a  combination  of 
the  norms  of  both  transfer  functions  Such  a  controller  attempts  to  optimize 
both  the  sensitivity  and  robustness  of  the  closed  loop  system.  The  relative 
importance  of  the  sensitivity  and  robustness  objectives  are  controlled  in  a 
frequency  dependent  fashion  by  multiplying  the  two  transfer  functions  by 
H°°  weighting  functions,  W\(z)  and  W2(z).  A  good  description  of  the  H°° 
design  methodology  may  be  found  in  [Fra87]  [DFT92],  and  details  of  its 
particular  application  to  the  pitch  control  problem  are  given  in  [EOT92]. 

In  short  the  design  specification  is  translated  into  the  requirement  that 
the  controller  should  be  chosen  to  stabilize  the  plant  and  minimize  the  value 
of  the  norm 


(  W2G{l  +  FG)~l 
Wi(l  +  FG)~l 


(5) 


where  W\  and  W2  are  the  weighting  functions  referred  to  earlier  in  the 
previous  paragraph.  If  the  Youla  parameterization  (4)  is  substituted  for  G 
in  (5)  then  the  objective  is  transformed  into  the  norm  of  the  affine  expression 


(  F2(z)(Xl(z)+F2(z)Q(z))W2(z)  \ 

V  F2(z)(X2(Z)  -  Fx{z)Q{z))Wi{z)  j  J  l  J 

The  design  problem  is  to  find  a  function  Q  in  H00  that  produces  a  value 
for  (6)  within  a  predefined  upper  bound.  This  problem,  which  is  related 
to  the  classical  Nehari  problem  has  been  at  the  center  of  frequency  domain 
approaches  to  the  H00  control  system  design  since  the  early  work  of  Francis 
et  al.  [FHZ84],  [FZ84],  The  authors  of  [EOT92]  and  [FY94]  use  operator 
theoretic  techniques  to  tackle  the  Nehari  Problem,  the  method  that  this 
paper  presents  starts  from  a  different  point  of  view.  If  the  affine  expression 
in  (6)  is  labeled  P,  then 


PM  -  (  F2(z)Xl(z)W2(z) 
[  ’  V  F2(z)X2(z)Wx(z) 


+  Q(z ) 


F2(z)2W2(z)  \ 

-F2{z)Fi{z)W1(z)  )' 


which  may  be  rewritten  as, 

(  F2(z)Xx(z)W2(z)  \ 
V  F2{z)X2{z)Wx{z)  ) 


1  P(z)  +  Q(z) 


—F2(z)2W2(z)  ) 
F2(z)Fi(z)Wi(z)  )' 


(7) 
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and  the  design  problem  may  be  viewed  as  the  problem  of  solving  a  linear 
Diophantine  equation  with  the  added  requirement  that  a  norm  on  the  so¬ 
lution  should  be  minimized.  More  specifically,  one  must  find  Q,  an  H°° 
function,  and  P,  a  matrix  of  H00  functions  which  together  solve  (7)  with 
the  added  constraint  that  the  L00  norm  for  the  largest  singular  value  of  P 
should  lie  within  a  specified  bound.  An  optimal  solution  is  one  for  which 
the  norm  on  P  is  a  minimum  with  respect  to  the  values  that  this  norm  takes 
over  all  solutions  P  to  Equation  (7).  The  algorithm  presented  here  will  not, 
in  general,  achieve  that  minimum,  and  we  demonstrate  in  Section  4  that  the 
discrepancy  can  be  large.  Equation  (7)  resembles  the  Bezout  equation  (3) 
which  is  also  a  linear  Diophantine  equation,  and  the  same  basic  technique 
is  used  to  compute  solutions  for  both  of  them. 
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3  Solving  Linear  Diophantine  Equations 

Let  /i,  /2  and  h  be  three  functions  in  H°°.  This  section  addresses  the 
problem  of  finding  solutions  hi  and  h2  in  H°°  for  the  equation 

fi9l  +  /2S2  =  h.  (8) 

Equation  (8)  subsumes  the  Bezout  equation  (3),  and  if  the  objective  function 
(6)  were  replaced  by  a  scalar  objective  function,  it  would  subsume  the  linear 
equation  (7)  that  arises  in  the  Nehari  problem  as  well.  If  f\  and  /2  are  outer 
functions,  which  is  to  say  that  they  possess  multiplicative  inverses  in  H°°, 
then  the  solution  is  easy.  A  family  of  solutions  with  parameter  g  an  H°° 
function  is  formed  by  setting  g\  =  gh /{" 1  and  g<i  =  (1  —  When  f\ 

and  /2  have  zeros  in  the  right  half  plane,  the  inverses  no  longer  exist  and 
this  method  breaks  down,  but  if  the  requirement  that  the  solutions  be  in 
H00  is  temporarily  relaxed  then  bounded  (but  not  analytic)  solutions  may 
be  found  as  follows. 

Let  0  be  a  bounded  function  on  the  half  plane  with  the  property  that  the 
zero  set  of  /1  is  bounded  away  from  the  support  of  </,  and  the  the  zero  set  of 
/2  is  bounded  away  from  the  support  of  1  —  0.  For  such  a  function  to  exist 
some  restriction  needs  to  be  placed  on  the  functions  f\  and  /2,  for  instance, 
it  is  necessary  that  they  should  have  no  common  zeros.  Bounded  solutions 
to  Equation  8  can  be  constructed  in  a  piecewise  fashion  by  taking  g\  —  0 
outside  the  support  of  <p,  and  g\  =  on  the  support  of  4>  and  for  the 

second  function  <?2  =  0  outside  the  support  of  1  -(f)  and  <?2  =  (1  —  0)/i/2'1 
on  the  support  of  1  —  cf>. 

Observe  that  if  e  is  a  bounded  function  on  the  right  half  plane,  then  the 
two  functions  — e f\  and  e /2  satisfy  the  relation 

{-e/i>/2  +  (e/2)/i  =  0. 

With  this  in  mind  the  step  from  the  bounded  solutions  g\  and  j2  to  H°° 
solutions  can  be  made  if  a  bounded  function  e  can  be  found  such  that  the 
functions 


91  =  9i  +  e/2  (9) 

92  =  92  ~  e/i 

are  are  both  in  H°°.  In  this  section  it  is  shown  that  if  the  function  (j>  is 
chosen  appropriately,  a  suitable  function  e  can  be  calculated  as  a  solution 
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to  a  first  order  partial  differential  equation.  The  H00  norms  of  the  solutions, 
gi  and  <72,  depend  on  the  choice  of  <f>  and  e.  In  the  case  of  the  Nehari 
Problem  described  in  section  2,  a  bound  on  the  norm  ||P||oo  is  obtained 
by  separately  bounding  quantities  that  correspond  to  g\  and  e  in  equation 
(9).  As  a  consequence,  success  in  calculating  solutions  with  small  ||P||oo 
norms  depends  critically  on  controlling  the  norm  of  the  function  e.  Section 
4,  which  presents  a  solution  to  the  problem  posed  in  Section  2,  contains 
further  discussion  of  this  point. 

3.1  Algebraic  Reformulation 

The  appropriate  setting  in  which  to  make  the  introductory  paragraphs  of 
this  section  precise  is  the  setting  of  homological  algebra.  This  setting  was 
first  presented  by  Hormander  [Hor67]  in  conjunction  with  the  corona  prob¬ 
lem,  and  has  been  used  by  Berenstein,  Taylor,  Struppa  and  Yger  [BS86], 
[BT80],  [Str83],  [BY89]  for  the  investigation  of  linear  Diophantine  equations 
in  algebras  over  rings  of  analytic  functions  of  bounded  growth.  In  the  sim¬ 
ple  cases  of  equations  that  arise  from  single-input  single-output  systems  the 
algebraic  formalism  reduces  precisely  to  the  equations  of  the  opening  para¬ 
graph.  A  good  introduction  to  the  analysis  used  in  this  chapter  is  [BG91], 
and  the  survey  [BS93]  describes  recent  applications  of  the  methods  that  are 
used  in  this  paper  to  other  problems  in  the  analysis  of  linear  operators. 

Let  R  denote  a  ring  of  functions  (or  distributions)  on  the  half  plane 
7i.  For  any  positive  integer  m  let  A (R)  denote  the  graded  module  over  R 
that  consists  of  functions  on  €  that  take  values  in  the  exterior  algebra  of 
antisymmetric  forms  on  an  m-dimensional  vector  space.  Let  A k(R)  denote 
the  homogeneous  elements  in  A(jR)  of  order  k.  A (R)  is  a  finite  free  module 

over  R ,  and  a  basis  element  will  be  denoted  by  . ik  with  the  indices 

ordered  1  <  i\  <  . . .  <  ik  <  m.  If  j  G  {ii, ... ,  ik)  then  the  symbol  e. 
denotes  the  basis  element  of  Afc_1  that  is  formed  by  deleting  the  index  j  from 
eu  if  j  £  {*1,  •  •  •  ,ik)  then  e.  7  .  =  0.  Suppose  that  /  =  {/1,  —  fm) 
is  a  finite  subset  of  R ,  then  the  operator  Pf  defined  on  A (R)  acts  on  a  basis 
element  as  follows 

k 

3= 1 

The  operator  Pf  forms  an  exact  sequence  over  the  homogeneous  sub- 
modules  A.k  called  the  Koszul  complex.  In  the  case  of  m  =  2  the  Koszul 
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complex  is  represented  by  the  diagram 


0 — >A2(R)-^Ua1(R)^Ua°(R) — >-0, 


(10) 


and  when  R  =  H00,  this  corresponds  exactly  to  the  concrete  algebraic  setting 
described  in  the  opening  paragraphs  of  this  section.  Finding  a  solution 
to  Equation  (8)  is  equivalent  to  inverting  the  operator  Pj  :  A1  (H°°)  —>■ 
A°(h°°)_  The  approach  to  inverting  the  operator  Pf  that  was  outlined  in 
the  introduction  was  first  to  find  an  inverse  in  a  larger  space  than  A1(H°°), 
and  then  add  an  element  from  the  image  of  Pf  :  A2(H°°)  ->  A1(H°°)  that 
will  return  the  solution  to  H00.  The  appropriate  ring  in  which  to  invert 
Pf  is  the  ring  of  distributions  which  will  be  introduced  with  the  following 
definitions  from  Hormander  [Hor67]. 

Denote  the  open  right  half  plane  by  PL,  its  closure  by  PL,  and  its  boundary, 
the  imaginary  axis,  by  dPL.  Let  fi  be  a  measure  with  support  on  PL,  then  a 
distribution  u  on  PL  satisfies  the  equation3 


du 

al 


=  P 


(11) 


if  for  any  continuously  differentiable  test  function  tp  with  support  compactly 
contained  in  PL, 


The  measure  dxdy  in  the  first  integral  is  the  Lebesgue  measure  on  C.  A 
distribution  u  that  satisfies  (12)  is  said  to  have  boundary  value  <p,  an  L°° 
function  on  the  imaginary  axis,  if  there  exists  U,  an  extension  of  u  to  PL, 
that  satisfies: 

nrr 

——  =  p  —  (pdz/2i.  (13) 

az 

Each  side  of  this  formula  is  to  be  interpreted  as  a  distribution  acting  on  test 
functions  supported  in  the  closed  half  plane  Pi  and  the  measure  <pdz/2i  is  a 
measure  on  €  with  support  on  the  imaginary  axis.  The  motivation  for  this 
definition  comes  from  Stokes’  Theorem. 

3  When  <D  is  considered  as  a  homeomorphic  to  R2  in  the  usual  way  (x,  y)  — >  x  +  iy  then 
the  operator  d/dz  is  expressed  in  local  real  coordinates  a s  d/dz  =  l/2(d/dx  +  id/dy ) 
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A  measure  p  in  fit  is  called  a  Carleson  measure  [Gar81]  with  Carleson 
constant  C  if 

n(S)  <  C  1{S)  (14) 

for  every  square  S  CfiL  with  a  side  of  length  l(S)  lying  on  an  interval  on  the 
imaginary  axis.  The  space  of  Carleson  measures  is  denoted  by  the  symbol 
C. 

Let  B  denote  the  ring  of  distributions  over  fit  with  boundary  value  in 
L°° ,  and  with  each  element  b  €  B  satisfying 


db 

dz 


=  A* 


(15) 


for  some  Carleson  measure  /x  in  fit.  This  is  the  ring  in  which  the  operator  Pf 
will  be  inverted.  The  differential  operator  djdz  :  B  ->  C  and  the  canonical 
injection  i  :  H°°  -»  B  form  an  exact  sequence 

0— tH°°-4£^4C— tO  (16) 


which  when  combined  with  the  Koszul  complex  (10)  gives  a  double  complex 
which  for  the  case  m  =  2  is  represented  by  the  commutative  diagram  of 
Figure  2.  For  notational  convenience  the  modules  A k(B)  and  A k(C)  are 
denoted  by  A*  and  A*  respectively.  Readers  familiar  with  complex  manifold 
theory  should  recognize  Equation  (16)  as  a  3  co-homology  sequence.  Since 
the  results  of  this  paper  are  restricted  to  analytic  functions  defined  on  the 
complex  half-plane,  it  is  not  difficult  to  avoid  introducing  the  language  of  the 
co-homology  of  differential  forms  —  the  appropriate  setting  for  analogous 
results  about  analytic  functions  of  several  complex  variables. 

The  next  theorem,  which  comes  from  [Hor67] ,  explains  how  the  complex 
is  used  to  provide  solutions  to  the  Diophantine  equations.  The  construction 
of  the  solutions  is  given  in  the  proof  which  is  repeated  here  for  the  sake  of 
completeness. 

Theorem  1  [Hormander] 

Suppose  that  the  following  conditions  are  satisfied: 

(i)  Let  s  take  the  values  0  and  1,  and  r  take  the  values  1  and  2.  If  h  €  A£ 
and  Pfh  =  0  then  the  equation  Pjg  =  h  has  a  solution  g  €  A*+1  with 
dg/dz  6  A^(  when  dh/dz  —  0. 

(ii)  dg/dz  =  p  has  a  solution  g  €  Af  for  every  p  €  A|. 
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0 


0 


0 


A2(H°°) 


A? 


d/dz 


Ai 


pf. 


0 

0 

0 


Figure  2:  double  complex  for  m  =  2 


Then  for  every  h  £  A3  with  dh/dz  =  0  one  can  find  g  £  A\  so  that  dg/dz  =  0 
and  Py<7  =  /i. 

Proof: 

The  result  follows  when  premises  (i)  and  (ii)  are  used  to  traverse  the  diagram 
in  Figure  2  as  follows. 

Suppose  that  h  £  A?  is  a  holomorphic  function  with  boundary  value  in 
L00;  that  is,  dh/dz  =  0  on  TL,  and  there  exists  a  function  H(y)  £  L°°(1R) 
such  that  for  almost  all  y  £  l R,  H(y)  =  lirn^jy  h(z)  when  the  limit  is  non- 
tangential  to  the  boundary.  Then  by  the  first  premise  there  exists  gl  £ 
Aj;  such  that  Pfg1  =  h  and  dgl /dz  £  A\-  Commutativity  implies  that 
Pfdg1  jdz  =  dj dz  Pjg 1  =  0,  so  again  by  the  first  premise  there  exists  g2  £  A 2 
such  that  Pfg2  =  dg1  jdz.  By  the  second  premise  the  equation  dg3 /dz  =  g2 
has  a  solution  g3  £  A2.  Let  g  =  gl  —  P/g3,  then 


% 

dz 


dg 1  p  dg 3 
dz  f  dz 
0 


and  Pfg  =  Pf(gl  -  Pfg3)  =  Pf{gx)  -  PfPfg3  -  P/g 1  =  h  as  required. 


□ 
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Before  Theorem  1  can  be  used  to  construct  solutions  to  Diophantine 
equations,  explicit  inversion  formulas  for  the  operators  P;  and  d/dz  satisfy¬ 
ing  premises  (i)  and  (ii)  need  to  be  presented.  It  was  this  need  for  inversion 
formulas  which  governed  the  definitions  of  the  spaces  B  and  C.  Formulae 
for  the  inversion  of  the  operator  Pj  that  are  based  on  the  work  in  [Hor67] 
with  only  small  modification  are  presented  next.  Only  the  case  of  m  =  2, 
the  case  needed  for  single-input  single-output  systems,  is  considered  here;  a 
more  general  situation  that  will  be  used  for  multi-input  multi-output  sys¬ 
tems  is  dealt  with  in  Section  5.  A  constructive  scheme  for  inverting  the 
Cauchy  Riemann  operator  d/dz  that  comes  from  the  more  recent  work  of 
Jones  [Jon80]  [Jon83]  is  presented  in  the  next  section. 

The  inversion  of  the  operator  Pf  :  A}  ->  is  dealt  with  first.  It 
turns  out  that  the  requirement  that  is  hardest  to  satisfy  is  the  requirement 
that  the  anti-holomorphic  derivative  of  the  inverse  should  be  a  bounded 
Carleson  measure.  To  overcome  this  problem  (which  in  fact  presents  a  major 
obstacle  in  the  proof  of  the  corona  theorem)  the  construction  is  based  on  an 
application  of  the  following  Lemma  from  [Hor67] 

Lemma  2  [Carleson  -  Hormander]  Let  fj  £  ft,  j  =  1, . . .  ,n,  and  assume 
that  for  some  c  >  0 

|/i(*)|  +  ...  +  |/»(*)|  >c.  (17) 

Then  for  sufficiently  small  e  >  0  one  can  find  a  partition  of  unity  <f>j  sub¬ 
ordinate  to  the  covering  of  ft  by  open  sets  ftj  =  {z  :  \fj(z)\  >  e}  such  that 
d(f>j/dz,  defined  in  the  sense  of  distribution  theory,  is  a  Carleson  measure  for 
all  j. 

This  Lemma  is  a  restatement  of  a  result  of  Carleson’s  original  paper  [Car62] 
in  which  he  directly  constructs  the  measure.  A  more  recent  account  of  the 
construction  is  given  in  Garnett’s  book  [Gar81].  The  difficult  part  of  the 
lemma  is  the  construction  of  a  partition  of  the  plane  into  two  sets  each  of 
which  contains  the  regions  of  the  plane  where  one  of  the  two  functions  f\  or 
fl  becomes  very  small.  In  general,  the  boundary  between  the  two  sets  will 
be  a  complicated  curve,  however  in  practice,  the  functions  f\  and  fi  may 
possess  some  regularity  that  allows  a  boundary  curve  to  be  easily  chosen. 
For  example,  when  solving  the  Bezout  equation  that  arises  from  the  example 
presented  in  Section  2,  f\  and  fi  are  the  functions  F\  (z)  =  e~TZ / (1  +  z)  and 
F2(z)  =  (1  -  z)/(l  +  z)  so  the  only  restriction  on  the  partition  is  that  it 
separate  the  point  z  =  1  where  ^(z)  =  0  from  the  regions  of  the  plane 
where  \z\  is  large  and  F\(z)  tends  to  zero,  and  a  simple  geometry  suffices. 
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The  partition  of  unity  from  Lemma  2  is  used  to  construct  a  left  inverse 
for  Pf  on  as  follows:  for  h  £  A®,  let 


4>i 

9i  =  h^~.  (18) 

Ji 

then  P/g  =  figx  +  fog^  =  h.  The  right  inverse  g  also  satisfies  the  premise  of 
Theorem  1,  for  if  dh/dz  —  0,  then  dg/dz  =  hf~1d(pj/dz  which  by  Lemma  2 
is  a  Carleson  measure. 

A  second  inversion  formula  is  needed  to  invert  the  operator  Pf  :  A|  — >  Ag. 
In  fact,  for  the  purposes  of  Theorem  1  it  suffices  to  invert  Pf  on  the  subspace 
consisting  of  measures  dgl/dz  where  g1  is  a  solution  of  Pfg1  =  h  for  some 
holomorphic  function  h.  In  this  case  the  problem  may  be  written  down 
explicitly  as  a  redundant  set  of  equations  for  the  coefficient  of  g2 

312/2  =  h/fi  d<j>i/dz 
-312/1  =  h!h  dfa/dz 


and  since  /2  and  /1  are  both  holomorphic  functions  with  magnitude  bounded 
away  from  zero  on  the  support  of  d(f>i/dz  a  solution  is  given  by 


d<t>\/dz. 


3.2  Constructing  Bounded  Solutions  to  the  Inhomogeneous 
Cauchy  Riemann  Equation 

In  the  preceding  section  the  construction  of  solutions  to  the  Bezout  Equation 
has  been  reduced  to  two  steps:  the  construction  of  the  partition  of  unity  <f)j, 
and  the  construction  of  bounded  solutions  to  the  Cauchy  Riemann  equation. 
This  section  describes  a  technique  devised  by  P.  Jones  [Jon80]  for  solving  the 
Cauchy  Riemann  equation;  the  presentation  is  based  on  the  account  given 
in  Garnett  [Gar81]. 

The  problem  that  needs  to  be  solved  is:  given  /i,  a  Carleson  measure 
on  the  right  half  plane,  find  a  distribution  b  with  bounded  boundary  values 
that  satisfies 

db/dz  =  g. 

The  solution,  which  is  based  on  a  Green’s  function  argument,  has  three 
stages:  the  measure  g  is  approximated  by  a  sequence  of  measures  gj  which 
converge  weakly  to  g,  each  gj  being  supported  on  a  finite  set  of  points;  the 
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support  of  each  measure  ytj  is  partitioned  in  such  a  way  that  the  pseudo- 
hyperbolic  distance4  between  any  two  points  in  the  same  partition  is  bounded 
from  below,  and  the  measure  fij  is  subdivided  into  a  corresponding  sum 
Yj  A tj  each  ^  having  support  on  a  distinct  set  in  the  partition;  finally  the 
Cauchy  Riemann  equation  is  boundedly  solved  for  each  yJj  and  these  solu¬ 
tions  are  summed  to  form  the  approximate  solution  bj.  The  whole  procedure 
is  performed  in  such  a  way  that  the  sequence  of  solutions  bj  is  a  uniformly 
bounded  sequence  of  functions  in  H00. 

Before  the  solution  is  discussed  in  detail  the  fundamental  solution  to  the 
Cauchy  Riemann  operator  d/dz  is  introduced,  and  a  result  about  interpo¬ 
lating  Blaschke  products  is  recounted.  Let  D  C  C  be  an  open  domain  with 
C1  boundary  that  contains  the  origin  z  =  0.  The  fundamental  solution  to 
the  operator  d/dz  on  D  is  a  distribution  b  that  satisfies  the  identity 

-  [  b(z)-^~-dxdy  =  <f>( 0) 

J  D  OZ 

for  any  C00  function  <f>  with  support  compactly  contained  in  D  [Hor90]. 
In  this  formula  the  integral  on  the  left  hand  side  of  the  identity  should 
be  interpreted  as  the  action  of  the  distribution  on  a  test  function.  The 
fundamental  solution  is  computed  as  follows.  Suppose  that  (j)  is  an  arbitrary 
C°°  function  with  support  compactly  contained  in  D.  Let  U  C  D  have  Cl 
boundary  and  contain  the  support  of  (f>  in  its  interior.  Consider  the  function 
</>(C)/C,  Stokes’  theorem  gives 


Because  cj)(z)  =  0  on  the  boundary  of  U,  the  first  boundary  integral  is  zero, 
and  as  e  -¥  0  the  second  integral  approaches  the  limit  i2ir<p(0).  Consequently 
a  fundamental  solution  for  d/dz  is  given  by  the  distribution  h(z)  =  1/(-kz). 

If  z  =  x  +  iy  is  a  complex  number,  then  the  real  conjugate  of  z  is 
defined  to  be  the  number  z  =  —x  +  iy.  The  need  for  this  usage  results  from 
considering  Laplace  transforms  of  system  operators;  the  Laplace  transform 
of  a  bounded  causal  system  gives  a  transfer  function  which  is  analytic  in 
the  right  half  plane,  so  in  places  where  a  complex  conjugate  z  occurs  in  the 

4The  pseudo-hyperbolic  distance  between  two  points  in  the  half-plane  is  defined  as 
p(zi,Z2)  =  \zi  -  Z2\l\z\  -Z2 1 
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analysis  of  functions  analytic  in  the  upper  half  plane,  it  will  be  natural  to 
substitute  the  real  conjugate  z.  For  instance,  given  a  set  {(j  =  +irjj,  > 
0}  that  satisfies  the  condition 


E 


1  +  ICnl2 


<  00, 


a  Blaschke  product  with  zeros  (j  is  defined  by  the  expression 


Kj  jj.£ _ 0 

0-i  2  -  o 


The  factors  |0  —  l|/(0  —  1)  ensure  that  the  product  converges  when  the 
sequence  {|0'|}  is  unbounded,  and  for  finite  zero  sets  they  may  be  omitted. 
A  finite  number  of  zeros  at  C  =  1  may  be  introduced  into  the  Blaschke 
product  separately  in  the  factor  in  front  of  the  product  sign. 

Let  Bo (z)  be  a  Blaschke  product  with  a  zero  set  (0  =  0  +  ir)j}  that 
satisfies  the  condition 


n 


Ck _ 0 

Ck  —  0 


>  s  >  o, 


(19) 


then  the  inverse  1/B0(z)  is  an  analytic  function  except  on  the  zero  set  (0} 
and  is  given  by  the  expression 


1/B0(*)  =  l  +  £ 

j 


i/gp(Cj) 

z-Cj 


If  l/Bo(z)  is  considered  as  a  distribution  on  B,  then  it  follows  from  the 
discussion  of  the  fundamental  solution  to  the  8  operator  that 

d 


dz 


(1/Bo) 


=  E 

i 


TV 


B'oiCj)  Cj 

~  52  > 


(20) 


where  1  <  \@j\  <  1/d. 

The  following  theorem  is  quoted  from  Jones  [Jon83]  Theorem  6. 
Theorem  3  [Jones] 

Suppose  {z/i}  is  a  sequence  of  points  in  the  half  plane  that  satisfies 


f  n 

J  k,kft 


Zk  ~  Zj 


zk  Zj 


>6>  0,  k  =  1,  2, . 
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Let  B(z)  be  the  the  Blaschke  product  with  zeros  at  the  points  { Zk }  and  Bj(z ) 
be  the  Blaschke  product  with  zeros  in  the  set  formed  by  removing  the  point 
Zj  from  the  set  {zk}-  Let  Ej(z)  be  the  function 


Ei(z)  =  CjBj(z) 


y  Uk  1 


where 


cj  —  —^{Bj(zj))  1  exp 


log  2/5 


E 

Vk<Vj 


Vk 


Zj  ~  zk 


Then  Ej(zk)  =  6jtk  and 


'Z,\Ej(z)\<(C0/6)  log(2 /S) 

3 


for  all  z  E  hi. 


(21) 


(22) 


This  theorem  is  an  instance  of  Carleson’s  interpolation  theorem  that  explic¬ 
itly  gives  the  form  of  the  interpolating  function.  The  bound  on  the  norm  of 
the  interpolating  function,  ( Cq/8 )  log(2/5)  is  optimal  in  5  up  to  the  multi¬ 
plicative  factor  Co- 

The  application  of  Jones’  interpolation  formula  requires  the  following 
lemma  which  is  extracted  from  the  proof  of  Carleson’s  interpolation  theorem 
in  Chapter  7  of  [Gar81]. 


Lemma  4  (Garnett  [Gar81]) 

Let  { Zj }  be  a  sequence  in  the  right  half  plane,  with  points  Zj  well  separated 
in  the  hyperbolic  metric,  i.e. 


- p(zk,Zj ) 


Zk  Zj 
Zk  -  Zj 


>  a  >  0, 


j  7^  k, 


and  suppose  that  there  exists  a  constant  A  such  that  for  every  square  Q  — 
{yo  <  y  <  yo  +  KQ),  0  <  x  <  l(Q)}, 


then 


w  n 


Zk  ~  Zj 
Zk  -  Zj 


>  5  >  exp  4(L4  ^1  +  2  log  — 
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The  bound  that  is  given  for  5  in  the  lemma  depends  on  the  points  {zk}  having 
a  minimum  spacing  a  in  the  hyperbolic  metric  and  on  the  measure 
being  a  Carleson  measure  with  Carleson  constant  A.  Unfortunately  the 
generality  of  the  theorem  means  that  the  bound  derived  will  be  conservative 
for  many  specific  examples.  This  is  particularly  true  of  examples  such  as  the 
one  presented  in  this  paper  in  which  the  measures  have  easily  recognizable 
additional  structure.  Additional  information  about  the  distribution  of  the 
points  {z^  could  well  be  used  to  derive  a  less  conservative  estimate. 

The  next  two  lemmas  contain  the  constructive  solution  to  the  Cauchy 
Riemann  equation  that  is  presented  in  Chapter  8  of  [Gar81].  The  proofs 
closely  follow  the  work  cited,  but  are  given  here  because  they  contain  the 
algorithms  that  are  used  to  compute  actual  solutions.  Jones’  interpolation 
theorem  and  the  discussion  preceding  it  on  fundamental  solutions  provide 
the  basis  for  calculating  solutions  to  the  Cauchy  Riemann  equation  in  the 
following  simple  case. 

Lemma  5  (Garnett  [Gar81]) 

Let  Zj  be  a  finite  set  of  points  satisfying  (19)  and  let  p  =  ctjXj5Zj  with 
\afi  <  1.  Then  the  function 

b{z)  =  E(z)/B\(z)  (23) 

satisfies  db/dz  =  p  where  B\(z)  is  a  Blaschke  product  with  zeros  Zj,  and 
E(z)  is  a  function  that  is  analytic  on  the  right  half  plane  and  has  a  bound 
that  depends  on  the  choice  of  p  only  through  the  S  of  Equation  (19). 

Proof: 

Equation  (20)  states  that 

^BZz )  = 

and  that  the  coefficients  j3j  lie  within  the  uniform  bounds  1  <  \(3j\  <  1/5. 
An  application  of  Jones’  interpolation  theorem  produces  a  function 

E(z)  =  YJ<*ilPiEj{z)  (24) 

that  is  analytic  in  the  right  half  plane,  interpolates  the  values  oy / (3j  at  the 
points  Zj ,  and  is  bounded  on  the  imaginary  axis  by 

\E(z)\  <  (Co/S)  log(2 /S) 

in  which  Cq  is  an  absolute  constant. 
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The  result  follows  by  taking  b(z)  —  E{z) / B\(z) 


□ 

The  case  of  a  general  Carleson  measure  y,  is  tackled  by  constructing 
a  sequence  of  approximating  measures  {yn}  that  converges  (weakly)  to  y\ 
each  measure  in  the  sequence  is  supported  on  a  finite  set  of  points  and  has 
the  form  yn  =  YJ  ajXjSZj .  A  sequence  {&„}  of  uniformly  bounded  solutions 
to  the  equations  dbn/dz  —  yn  is  calculated,  and  the  restrictions  to  the 
imaginary  axis  {bn(iy)}  form  a  uniformly  bounded  sequence  of  L°°  functions. 
The  mapping 

/OC 

f{y)K{iy)  dy,  f  eLl 

-00 

associates  the  set  {bn(iy)}  with  a  set  of  uniformly  bounded  functionals  on  Ll, 
and  weak  compactness  of  the  unit  ball  ensures  that  {bn(iy)}  contains  a  sub¬ 
sequence  that  converges  in  the  weak*  sense  to  an  H°°  function  b(iy).  Choose 
one  such  subsequence,  and  relabel  it  {bn(iy)},  then,  if  ip  is  any  Co°(C)  test 
function  with  support  that  intersects  the  half  plane,  the  restriction  to  the 
imaginary  axis  ip(iy)  is  an  Lx  function,  and 

/oo  roo 

ip{iy)b{iy)dy  =  lim  /  ip(iy)bn{iy)  dy 
-OO  J- oo 

=  lim  t  [  i pbn(z)  dz 
n-ioo  i  Jdn 

=  lim  —  /  d{ipbn)  A  dz 

n-*00  1  Jn 

1  f  dip  _  f  dbn 

=  lim  -  /  bn—dz/\dz+  /  ip-r—dzf\dz 
n-yoo  i  Jn  dz  Jn  OZ 

f  dip  f  dbn 

=  lim  2  /  bn~  dx  A  dy  +  /  ip-—rdxAdy 
n-*°°  Jn  dz  Jn  dz 

Rearranging  the  left-most  and  right-most  sides,  and  substituting  yn  =  dbn/dz 
gives  the  equation 

f  d'lb  /  d'tb  \ 

JH  b- s  dx  A  dy  =  V’  a?/  “  i’Mbl-iv)dy 

which  provides  a  consistent  definition  for  b,  as  a  distributional  solution  to 
db/dz  =  /x  on  J-L  with  boundary  value  b(t) 
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Lemma  5  is  not  quite  enough  to  provide  the  sequence  of  solutions  {6n}; 
the  difficulty  is  that  the  bound  in  Lemma  5  depends  on  the  parameter  5 
which,  through  Lemma  4,  is  related  to  the  spacing  (in  the  pseudo-hyperbolic 
metric)  of  the  points  in  the  supporting  set  { Zj },  and  if  a  general  Carleson 
measure  is  going  to  be  approximated  by  a  sequence  of  measures  with  finite 
point  support,  then  the  spacing  of  the  points  in  the  support  of  the  approxi¬ 
mating  measures  will  decrease  to  zero  as  the  approximations  converge.  What 
is  needed  is  a  method  for  decomposing  the  approximating  measures  in  such 
a  way  that  the  spacing  between  points  of  support  for  each  element  of  the 
decomposition  remains  large,  yet  the  sum  of  the  Carleson  constants  of  the 
elements  in  the  decomposition  remains  constant.  The  next  lemma  uses  this 
approach  to  produce  a  method  for  solving  the  equations  dbn/dz  =  pn  with 
a  uniform  bound  on  the  sequence  of  solutions  bn. 

Lemma  6  (Jones-Garnett  [Gar81]) 

Let  p  =  J2jLi  ajxj^zj  be  a  measure  supported  on  the  finite  set  {zj  =  Xj  + 
iyj},  with  masses  ajXj  at  the  points  Zj,  and  with  Carleson  constant  N(p)  < 
C.  Then  there  exist  an  integer  N,  functions  bp(z),  and  a  function 

2  N 

K* )  =  H  Mz) 

p=i 

such  that  each  bp(z)  is  a  function  of  the  type  produced  in  Lemma  5,  db/dz  = 
p  for  a  measure  p  that  is  arbitrarily  close  to  p,  and  |6(zt)  |  <  KC  for  t  €  1R 
and  K  a  constant  independent  of  p. 

Proof: 

First  it  is  shown  that  p  may  be  approximated  arbitrarily  closely  by  a  new 
measure  p  of  the  form  C/N  Y  Xj8Zj .  The  support  of  p  is  the  same  as  the 
support  of  p ,  but  each  point  mass  Xj6Zj  may  be  repeated  a  finite,  and  possi¬ 
bly  large  number  of  times  in  the  new  sum.  If  N  is  chosen  to  be  a  sufficiently 
large  positive  integer,  the  coefficients  ctj  in  the  finite  sum  p  may  be  uni¬ 
formly  approximated  to  arbitrary  accuracy  by  a.j  «  nj/NC  in  which  rij  are 
positive  integers  and  C  is  the  Carleson  constant  of  p.  If  each  term  in  the 
sum  Yj  cxjXj5Zj  is  expanded  as 

Q 

ajXjSZj  ss  —(xjSzj  +  ...  (nj  times) .. .  +  XjSzfi 
then  a  renumbering  of  the  terms  in  the  summation  gives  the  approximation 

A*  ~  A*  —  *^2  xjSzj 
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From  here  on  no  distinction  will  be  made  between  the  measure  fj,  and  the 
approximation  /2. 

In  the  second  part  of  the  proof  a  systematic  method  of  decomposing  the 
measure  /i  is  established.  The  point  masses  XjSZj  are  distributed  amongst  a 
finite  number  of  sets  in  such  a  way  that  the  distance  between  any  two  points 
in  the  same  set  is  large  in  the  hyperbolic  metric. 


Figure  3:  Dyadic  Subdivision  of  the  Half  Plane 


Choose  a  square  Qo  with  supp/u  C  Qo  C  %  and  a  side  of  length  l(Q0 ) 
lying  on  the  imaginary  axis.  This  square  may  be  subdivided  to  form  a  dyadic 
sequence  of  squares  of  uniform  hyperbolic  size  as  follows  (Figure  3  illustrates 
the  construction).  Let  Qi,  Q2  be  the  two  adjacent  squares  that  comprise  the 
left  half  of  the  square  Qo-  Each  have  sides  of  length  l(Q)/ 2,  and  each  have 
one  side  on  the  imaginary  axis;  continue  this  subdivision  process  inductively 
on  each  square  Qi  until  the  squares  Qin^  are  outside  the  support  of 

/i  for  some  n  (the  process  is  guaranteed  to  stop  because  the  support  of  /j,, 
which  is  finite,  is  compactly  contained  in  %).  Since  the  Carleson  constant 
of  n  is  fixed  to  the  constant  C,  a  simple  count  shows  that  the  right  hand 
section  of  any  dyadic  square  Q  can  contain  at  most  2 N  points  Zj.  This  allows 
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the  points  {zj}  to  be  partitioned  into  2 N  sets  {5n}  in  such  a  way  that  the 
spacing  between  any  two  points  in  the  same  set  is  uniformly  bounded  from 
below  by  a  =  1/3. 

The  sets  Sn  are  explicitly  defined  as  follows.  For  every  n,  let  Sn  =  {zj  : 
l(Qo)2~n~l  <  Xj  <  l(Qo)2~n}  and  order  the  elements  of  each  Sn  so  that 
Sn  =  { %k,n  +  Wk,n}  with 


Vk— l,n  —  Vk,n  —  yk+\,n- 


Then  the  set  {zj}  may  be  split  into  2N  sequences  hi .... ,  Yjyy  such  that 
the  points  in  each  Sn  are  evenly  distributed  between  the  Yr,  i.e.  if  Zj  = 
Xk,n  +  Wk,n  £  Sn  then  Zj  eYr  if  r  —  k  mod  2 N.  Now  suppose  that  P  C  T-L 
is  a  fixed  square  of  arbitrary  size  with  one  side  lying  on  the  imaginary  axis, 
let  Mn(P )  be  the  number  of  points  in  Sn  fi  P,  then  each  set  YT  fl  Sn  fl  P 
must  contain  fewer  than  1  +  Mn(P)/(2N)  points  Zj)  and 


E  (i+^)2-»(«3o) 

YrDP  n:SnC)P?<p  V  7 

OO  1 

n= 0  ZjgP 

<  4 l{P)+fi{P) 

<  5 1{P)  (25) 

Consider  the  sets  {Xp}  defined  by 

Xr  =  Yrn  U  Sn  (26) 

n  even 

and 

X2r+1  =  Yr  D  U  sn.  (27) 

n  odd 

then  the  measures  np  =  J2Zj €XP  xi$zj  satisfy  /i  =  C/N  p,p  and  up  to  the 
factor  C/N  provide  a  decomposition  of  /x  into  measures  with  well  spaced 
support  and  Carleson  constant  uniformly  bounded  by  5. 

A  bound  on  the  separation  between  points  of  support  is  arrived  at  by 
the  following  argument.  If  Zi  £  Sn  and  Zj  £  Sn- 2  then  p(zi,Zj)  >  1/3  by 
the  definition  of  the  sets  Sn.  On  the  other  hand,  if  Z{  and  Zj  are  in  the  same 
set  Sn  then  it  follows,  from  the  fact  that  the  top  half  of  any  Qj  contains 
at  most  2 N  points,  and  the  way  in  which  the  set  Yr  that  corresponds  to 
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Xp  was  constructed,  that  z*  and  Zj  must  be  separated  by  at  least  1  square 
of  length  l(Qo)2~n.  Consequently,  the  distance  between  Z{  and  Zj  must  be 
bounded  below  by 

2~n 

pizu  zj)  >  v/2-2n  +  8-2n 

>  1/3.  (28) 

Lemma  5  can  now  be  applied  to  the  measures  pp  to  produce  functions 
bp  that  satisfy 

dbp 

The  constant  S  in  Lemma  5  which  determines  the  bounds  on  the  norms  ||6p|| 
is  estimated  by  using  Lemma  4  and  the  inequalities  (25)  and  (28).  This 
gives  an  estimate  on  the  norms  ||6p|j  of 

IM  <  K, 

in  which  K  is  an  absolute  constant  that  is  independent  of  the  measure  pp. 
Let 

r  2  N 

H?)  =  Tfl Zhp(z)- 

P=  1 

Then  db/dz  =  C/N  X)  zj$Zj  =  P  and 

\\b\\  <  2CK,  (29) 

which  completes  the  proof. 

□ 

Lemma  6  provides  the  last  step  in  the  constructive  proof  of  the  following 
theorem. 


Theorem  7  (Garnett  [Gar81]) 

Let  n  be  a  Carleson  measure  with  Carleson  constant  N(/j.)  <  1.  Then  there 
is  a  distribution  b(z)  with  L00  boundary  value,  supported  on  H  such  that 

db 

¥z=^ 

and  the  boundary  value  satisfies  ||5||oo  <  C  for  C  a  positive  constant  in¬ 
dependent  of  the  choice  of  p.  Further,  there  exists  a  sequence  of  measures 
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{lik}  that  satisfy  the  criteria  of  lemma  6,  and  which  converge  weakly  to  p. 
The  corresponding  sequence  of  solutions  {bk(z)}  converge  in  a  distributional 
sense  on  TL  to  b[z),  and  have  boundary  values  that  converge  in  the  weak-star 
topology  to  b(iy)  on  the  imaginary  axis. 

The  reason  for  presenting  Lemma  6  in  such  detail  is  that  the  construction 
in  the  proof  provides  a  key  part  of  the  algorithm  that  is  used  to  compute 
solutions  to  the  Diophantine  equations  arising  from  the  control  problems.  In 
this  application  a  bound  on  the  norm  of  the  solution  to  the  equation  db/dz  = 
p  has  physical  significance,  and  consequently  a  tight  a  priori  estimate  of 
this  bound  would  be  valuable.  Unfortunately,  the  generality  of  the  methods 
presented  means  that  the  estimates  on  the  norms  that  can  be  obtained  from 
Lemmas  4,  5  and  6  are  too  conservative  to  be  of  practical  use. 


25 


4  Solution  to  Enns’  example 


The  theory  presented  in  the  previous  section  provides  a  practical  way  to  de¬ 
sign  linear  compensators  for  a  general  class  of  linear  time  invariant  systems. 
The  remainder  of  this  paper  provides  examples  that  illustrate  the  use  of  the 
theory.  This  section  continues  the  example  that  was  started  in  Section  2  by 
showing  how  a  compensator  is  calculated  for  the  plant  that  was  given  in  the 
example. 

The  plant  from  the  example  in  Section  2  is  described  by  the  transfer 
function  in  Equation  (1), 

e~TZ 
crz  —  1 


Values  are  ascribed  to  the  parameters  for  the  numerical  calculations:  er  takes 
the  value  a  =  1,  making  the  open  loop  system  unstable,  and  r  takes  the 
values  r  =  0.06  or  r  =  0.37,  the  first  is  a  small  delay  which  has  little  influence 
on  the  behavior  of  the  open  loop  system,  and  the  second  is  a  large  value 
for  the  delay  that  makes  the  problem  of  robust  stabilization  significantly 
more  challenging.  The  transfer  function  of  a  stabilizing  controller  is  given 
in  terms  of  an  H°°  parameter  Q  by  the  bilinear  function  in  Equation  (4). 
Evaluating  this  expression  requires  solutions  to  the  Bezout  equation  (3)  and 
the  Diophantine  equation  (7),  the  Bezout  equation  is  dealt  with  first. 


4.1  Computing  solutions  to  the  Bezout  equation 

Using  the  notation  of  Section  3.1  the  Bezout  equation  may  be  re-written 
as  Pfg  —  h  in  which  h  is  the  constant  function  h(z )  =  1,  /  is  the  pair 
of  co-prime  factors  fi{z)  =  N(z)  and  fa  =  M(z),  and  g,  the  solution,  is 
a  function  in  Al(R)  with  components  g\  =  X\  and  <72  —  X2  which  are 
H00  functions.  From  theorem  1  the  solution  may  be  written  as  g  =  gl  — 
Pfg 3  in  which  gl  satisfies  the  equations  Pfg1  =  h  and  dgl/dz  =  Pfg 2  for 
g2  a  Carleson  measure  on  C,  and  g 3  satisfies  dg3/dz  =  g 2  for  the  same 
measure  g2.  Substituting  the  notation  of  the  problem  gives  the  following  set 
of  equations: 

Xx{z)  =  Xx  (z)-b(z)F2(z)  (30) 

X2(z)  =  X2(z)  +  b{z)Fi{z) 

in  which  X\  and  X2  satisfy  the  equation 

Xl(z)Fl(z)+X2(z)F2(z)  =  l,  (31) 
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and  b(z )  =  g\2  =  -c/21  is  the  coefficient  function  of  the  2-form  g2  and  a 
distributional  solution  of 

db 
dz 

The  inhomogeneous  term  p  in  equation  32  is  a  Carleson  measure  that  is 
supported  on  the  half  plane  and  that  satisfies 


dXx 

dz 

dX2 

dz 


gF2 


—gF\. 


(33) 


Lemma  2  guarantees  the  existence  of  choices  for  Xx,  X2  and  p  that  satisfy 
equations  (31)  and  (33).  In  fact,  since  the  function  Fx(z)  =  e~hz / (z  +  1)  is 
bounded  away  from  zero  on  any  set  compactly  contained  in  the  right  half 
plane,  and  the  function  M(z )  =  (1  —  z)/(l  +  z)  has  a  single  zero  at  the 
point  z  =  1  and  is  bounded  away  from  zero  on  any  set  that  excludes  a 
neighborhood  of  that  point,  a  partition  of  unity  that  satisfies  the  condition 
in  lemma  2  is  the  following: 


<t>  1(2)  = 


1,  \z  —  1|  <  r 
0,  \z  —  lj  >  r 


M*) 


1)  |z  —  1|  >  r 
0,  |z  —  lj  <  r 


(34) 


With  this  choice  for  a  partition  the  solution  to  (31)  determined  by  equation 
(18)  is  Xx(z)  =  1/F\{z)  (f>i[z)  and  X2(z)  =  1/F2(z )  <p2{z).  Taking  anti- 
holomorphic  derivatives,  and  substituting  into  the  first  of  the  equations  (33) 
yields  the  measure 

H  =  l/(Fx(z)F2(z))d<f>l/dz.  (35) 

The  meaning  of  the  expression  d(j>\/dz  may  be  elucidated  by  mollifica¬ 
tion.  Let  12  denote  the  support  of  <j>\ ,  and  let  'ipk  be  a  sequence  of  positive 
C°°  functions  supported  on  connected  neighborhoods  of  the  origin,  and  with 
the  property  that  diameter (supp^fe)  -*  {0}  as  k  — >  00.  For  each  k  the  C°° 
function  cpk  =  4>  1  *  '’Pk  is  a  mollification  of  <p\\  for  sufficiently  large  k  it  is 
supported  on  a  region  slightly  larger  than  12,  and  takes  the  constant  value  1 
on  a  region  slightly  smaller  than  12.  Let  Dk  denote  the  support  of  the  C°° 
function  rjk  =  dcpk/dz.  It  follows  that  is  a  tubular  neighborhood  of  the 
boundary  of  the  support  of  and  that  (f>k  takes  the  value  1  on  the  interior 
part  of  dD  (the  boundary  of  D),  and  the  value  0  on  the  exterior  part  of  dD. 
The  function  %  may  be  interpreted  as  a  complex  valued  measure  on  C  in 
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the  following  sense.  Let  x  be  a  compactly  supported  C°°  function,  then  an 
application  of  Stokes  theorem  gives 


dxdy 


d(pk  , 

X~z=- dz  A  dz 
lDk  dz 


= 

2  JD_ 

if-  if-  dx 

=  o/  <t>kXdz-~  4>k—dzAdz 
2  JdDi.  2  Jz>.  dz 


(36) 


As  k  -¥  0  the  sequence  (f)k  converges  to  <f> i  in  the  topology  of  the  space  of  dis¬ 
tributions,  the  area  of  the  region  Dk  converges  to  0,  and  the  boundary  dDk 
converges  to  the  set  dCl.  Since  the  function  dx/dz  is  uniformly  bounded, 
the  second  integral  in  (36)  converges  to  0,  and  the  first  integral,  which  only 
has  a  contribution  from  the  interior  part  of  dDk ,  converges  to  —i/2  JdQ  x  dz. 
The  negative  sign  on  the  contour  integral  is  a  consequence  of  the  orientation 
of  the  boundary  dDk.  The  expression  d<f>\/dz  is  interpreted  as  a  measure 
supported  on  the  set  50,  which  acts  on  a  C°°  function  x  by 


The  shape  of  the  region  O  is  arbitrary  provided  that  the  partition  of 
unity  that  it  determines  satisfies  the  condition  in  lemma  2.  For  the  actual 
computation  of  the  solutions  to  the  Bezout  equation,  O  was  chosen  to  be 
a  circle  for  the  pragmatic  reasons  that  it  is  a  simple  curve  to  describe  and 
that  it  seems  to  give  reasonable  results.  The  radius  for  the  circle  was  chosen 
with  a  view  to  keeping  the  Carleson  constant  of  the  measure  /i  small.  It 
was  found  that  r  =  0.7  is  a  suitable  value  when  the  delay  in  the  plant  takes 
either  of  the  values  r  =  0.06  or  r  =  0.37. 

The  only  remaining  step  in  determining  solutions  to  Equation  (3)  is 
the  calculation  of  a  solution  to  the  Cauchy  Riemann  equation  (32).  The 
algorithm  presented  in  the  proof  of  Lemma  6  provides  a  way  to  calculate 
numerical  approximations  to  a  solution  of  the  equation.  The  approximations 
converge  in  the  weak-star  sense;  this  is  a  natural  topology  in  the  context 
of  system  transfer  functions.  In  accordance  with  the  proof  of  Lemma  6, 
the  measure  /i  is  approximated  by  a  finitely  supported  measure  constructed 
as  a  sum  of  point  masses  uniformly  distributed  on  the  support  of  //,  and 
the  algorithm  presented  in  Lemma  6  is  used  to  calculate  an  approximate 
solution  for  db/dz  =  fj,. 
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4.2  Computing  solutions  to  the  Nehari  equation 

The  Diophantine  equation  (7)  is  solved  using  the  same  method  as  was  used 
for  the  Bezout  equation  with  only  small  modifications.  The  similarity  in 
form  between  the  Nehari  problem  and  the  Bezout  equation  becomes  appar¬ 
ent  when  Equation  (7)  is  rewritten  as 


A  =  IP  +  BQ  (37) 

in  which 

_  (  f2xxw2  \  (  -f%w2  \ 

F2X2Wx  J  {  f2f1w1  j  ■ 

For  purposes  of  comparison  the  choice  of  weighting  functions  Wx  and  W2  is 
the  same  as  that  used  in  [EOT92]: 


Wi 

w2 


2  1  +  2 
1  +  10  z 

0.2  . 


Equation  (37)  may  be  rewritten  again  in  the  notation  of  Section  3.1  as  P/g  = 
h.  Now  h  =  A  is  a  vector  with  two  H°°  components,  /  has  components 
/i  =  1  and  f2  =  B,  a  vector,  and  the  solution  g  has  two  components, 
gi  =  P,  a  vector  with  two  H°°  components,  and  g2  =  Q  an  H°°  function. 
The  same  formalism  that  was  used  for  the  Bezout  equation  yields  the  set  of 
equations: 

P  =  P  —  b(z)B  (38) 

Q  =  Q  +  b{z) 

where  P  and  Q  satisfy 

P+QB=A 

and  b(z)  is  a  distributional  solution  of 

db 

where  g,  is  a  Caxleson  measure  that  satisfies  either  of  the  equivalent  equations 


dP 

dz 

dQ 

dz 


gB 

-g. 
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Figure  4:  Contour  for  measure  in  Nehari  Problem 


At  this  point  a  problem  becomes  apparent.  Consider  P  to  be  a  2  x  1 
matrix  with  H°°  entries,  then  the  optimal  solution  to  the  control  problem 
is  the  solution  that  minimizes  the  L°°  norm  of  the  largest  singular  value  of 
P(iy).  The  problem  is  that  the  algorithms  based  on  the  methods  presented 
in  section  3  can  only  guarantee  an  upper  bound  for  the  norm  on  P,  they  do 
not  give  information  about  how  close  the  computed  solution  is  to  an  opti¬ 
mal  solution.  The  upper  bound  on  the  solution  is  obtained  by  the  triangle 
inequality  from  bounds  on  the  size  of  the  two  terms  in  the  right  hand  size 
of  equation  (38).  The  bound  on  the  size  of  the  first  term  P  is  controlled 
directly  by  the  choice  of  the  distributional  solutions  P  and  Q,  and  the  bound 
on  the  second  term,  bB  is  controlled  by  the  the  Carleson  constant  for  the 
measure  y  and  the  a  priori  bound  that  theorem  7  imposes  on  b  the  solution 
to  the  Cauchy  Riemann  equations.  In  the  example  calculation  described  in 
the  following  paragraphs  ad  hoc  methods  are  used  in  the  selection  of  the 
distributions  P  and  Q.  The  results  indicate  that  better  methods  of  selec¬ 
tion  need  to  be  combined  with  improved  solutions  for  the  Cauchy  Riemann 
equations  if  the  goal  of  reliable  computation  of  near  optimal  controllers  is 
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to  be  achieved. 

The  distributions  P,  Q  are  chosen  to  minimize  the  L°°  norm  of  P  and 
the  Carleson  constant  for  /x.  The  equation  (^4  -  QB )  •  B  =  0  has  a  solution 
Q  =  A  ■  B/B  ■  B  which  is  a  bounded  analytic  function  away  from  the  zeros 
of  the  denominator  B  ■  B,  and  unbounded  in  a  neighborhood  of  these  zeros. 
Let  be  a  region  in  the  complex  plane  that  contains  a  neighborhood  of  the 
zeros  of  B  and  is  bounded  by  an  absolutely  continuous  closed  curve,  and  let 
4>i(z)  =  1  when  z  e  Cl  and  <j>\(z)  —  0  when  z  £  Cl,  and  let  Q  and  P  be  the 
bounded  distributions: 


A  •  B 
B-B 


<t>  i 


P 


A-QB 


A  — 


A-B 
B  ■  B 


B<f>x. 


(39) 

(40) 


Outside  the  region  Cl  the  measure  P  is  the  pointwise  orthogonal  projection 
of  A  onto  Bl.  Substituting  in  the  values  for  A  and  B  gives 


Q 

Pi 

K 


w2f2xx  -  W?FxX2  , 
W$F$  +  WlF'(  <Pl 


F2XxW2(1  -  <f> i)  + 
F2X2Wx(1  —  (f>\)  + 


F2FxWlW2  L 

w2F$  +  wxwrl 
f2wxw‘£ 

w£f£  +  wxw£  011 


and  for  the  measure  /i  acting  on  a  C°°  function 

i  f  w£f2xx-w?f2x2  , 
~  2  Jen  W£f£  +  W£f£  xdz 


(41) 


The  criteria  for  choosing  the  region  O  are  that  it  should  include  the  ze¬ 
ros  of  the  function  W£M2  +  W2N2  and  that  the  L°°  norm  of  P(iy )  and 
Carleson  constant  of  the  measure  fi  should  be  minimized.  For  the  re¬ 
sults  presented,  the  region  Cl  was  chosen  by  plotting  the  weighting  function 
(W2MY  —  WxNX)/(W2M2  +  WxN2)  and  choosing  by  inspection  a  contour 
that  includes  the  singularities  of  the  function,  yet  keeps  the  Carleson  con¬ 
stant  of  /x  small.  The  chosen  contour  is  illustrated  in  Figure  4;  the  contour 
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is  superimposed  on  a  logarithmically  scaled  contour  plot  of  the  magnitude 
of  the  weighting  function  in  Equation  (41).  Once  the  contour  is  chosen,  the 
Carleson  measure  in  (41)  can  be  computed  and  the  inhomogeneous  Cauchy 
Riemann  db/dz  =  fx  may  be  solved.  The  solution  b  is  substituted  into  (38) 
to  give  the  Youla  parameter  Q  and  the  transfer  function  matrix  P  that 
constitute  a  sub-optimal  solution  to  the  Nehari  Problem. 

4.3  The  Algorithm  and  its  Implementation 

This  section  contains  a  summary  of  the  algorithm  which  indicates  the  ex¬ 
act  order  of  the  steps  taken  in  obtaining  solutions,  and  where  the  various 
formulae  in  Sections  2-4  are  used. 

The  data  for  the  problem  are  the  transfer  function  F  given  in  equation 
(1),  and  the  weighting  functions  that  appear  in  the  norm  (5).  The  solution 
algorithm  can  be  divided  into  two  parts,  one  part  to  calculate  the  solutions 
X\  and  X2  to  the  Bezout  Identity  (3),  and  the  other  to  calculate  solutions 
to  the  Diophantine  equation  (7)  associated  with  the  Nehari  Problem. 

Parti:  The  Bezout  Identity 

1  Factor  the  transfer  function  —  Equation  (2). 

2  Choose  a  partition  of  H  that  separates  the  zeros  of  the  factors  F\  and 
F2  —  Equation  (34). 

3  Define  the  Carleson  measure  fx  for  the  Bezout  equation  (35). 

4  Compute  a  discrete  approximation  to  the  Carleson  measure  as  a  finite 
sequence  of  complex  valued  point  masses  at  points  on  the  support  of 
fx  —  Theorem  7. 

5  Group  the  points  according  to  the  algorithm  presented  in  the  proof  of 
Lemma  6. 

6  Use  Equation  (21)  of  Theorem  3  and  equations  (23)  and  (24)  of  Lemma 
(5)  to  compute  solutions  to  the  inhomogeneous  Cauchy  Riemann  equa¬ 
tion  (32). 

7  Use  Equation  (30)  to  calculate  the  solutions  to  the  Bezout  equation 
at  points  z  £  PL . 
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In  practice,  the  solution  is  calculated  over  a  given  set  of  points.  Steps  1 
-  5,  which  are  independent  of  the  point  z  at  which  the  solution  is  to  be 
calculated,  are  performed  first,  and  then  the  calculations  in  steps  6  and  7 
are  performed  at  each  point  2  in  the  given  set. 

Part2:  The  Diophantine  Equation  for  the  Nehari  Problem 

8  Choose  a  partition  of  T-L  so  that  the  Carleson  constant  for  the  measure 
in  Equation  (41)  is  minimized.  The  Carleson  constant  for  a  measure 
is  defined  in  Equation  (14). 

9  Define  the  Carleson  measure  /j,  for  the  Nehari  problem  by  Equation 
(41). 

10  Follow  steps  4-7  above  to  compute  the  solutions  to  the  inhomo¬ 
geneous  Cauchy  Riemann  equation  associated  with  the  measure  con¬ 
structed  in  step  9.  This  step  requires  the  evaluation  of  X\  and  X-i  at 
each  of  the  points  in  the  discretization  of  the  Carleson  measure  /i  from 
step  9. 

11  Use  Equation  (38)  to  evaluate  the  solutions  P  and  Q  on  points  z  G  H. 

The  whole  algorithm  for  computing  the  solutions  is  encoded  in  computer 
software.  The  heart  of  the  computation  is  performed  by  two  C  programs. 
The  first  program  computes  the  approximation  to  the  Carleson  measure,  and 
the  second  computes  the  uniformly  bounded  approximations  to  the  solution 
of  the  Cauchy  Riemann  equation.  Matlab  routines  that  call  the  C  programs 
calculate  the  transfer  functions  for  the  controller  on  a  set  of  points  on  the 
imaginary  axis,  and  produce  the  graphical  output. 

The  first  C  program  produces  a  discrete  approximation  to  the  measure  /i 
and  orders  the  points  of  support  according  to  the  decomposition  described 
in  Lemma  6.  The  process  occurs  in  a  number  of  discrete  steps.  The  first 
step  is  to  calculate  a  discrete  approximation  to  the  measure.  The  measure 
fj,  is  supported  on  a  curve  in  the  right  half  complex  plane.  In  the  example 
presented  here,  in  both  the  case  of  the  solution  to  the  Bezout  equation  and 
the  case  of  the  solution  to  the  Diophantine  Equation  associated  with  the 
Nehari  Problem,  this  support  is  compactly  contained  in  the  open  right  half 
plane,  and  the  weighting  function  is  fairly  regular.  Consequently  it  suffices 
to  choose  the  supporting  points  of  the  discretization  to  be  equally  spaced  on 
the  supporting  curve.  A  case  of  non-compact  support  would  require  that  the 
measure  be  restricted  first  to  a  compact  subset  of  the  plane,  and  the  added 
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complexity  of  non  uniform  spacings  in  the  discretization  may  be  justified 
in  some  instances  by  an  increase  in  the  efficiency  of  the  approximation. 
The  discretization  of  the  measure  is  stored  as  a  list  of  ordered  pairs  of 
complex  numbers,  the  first  element  in  each  pair  is  the  point  of  support 
Zk  —  %k  +  iyk ,  and  the  second  is  the  value  of  the  weight  assigned  to  each 
point  —i( Axk  +  iAyk)  /  (2F{ (z/t)F2(z^)).  Axk  +  iAyk  is  the  tangent  to  the 
supporting  curve  at  Zk  with  length  equal  to  the  arc-length  spacing  of  points 
in  the  discretization. 

The  second  part  of  the  measure  computation,  the  ordering  of  the  points, 
is  independent  of  the  particular  measure  being  approximated.  The  algo¬ 
rithm  proceeds  via  a  sequence  of  procedures  that  implement  the  steps  in  the 
proof  of  Lemma  6.  The  first  procedure  takes  the  list  of  weights  from  the 
discretization  and  and  produces  a  list  of  weights  of  uniform  size  by  splitting 
large  weights.  Given  a  small  constant  S  >  0  there  exists  an  integer  R  suffi¬ 
ciently  large  that  each  weight  may  be  written  a s  Wk  =  mkWk/R  where 
is  an  integer  and  1  <  |to/t|  <  1  +  <5.  So  the  weight  (zk,Wk)  is  split  to  form 
mk  weights  ( Zk,Wk/R )•  A  second  procedure  groups  the  weights  into  sets 
according  to  the  value  of  the  real  part  of  their  support,  the  weights  in  each 
band  are  sorted  by  the  imaginary  part  of  their  support,  and  these  sorted 
bands  are  passed  to  a  procedure  that  selects  weights  from  the  sorted  bands 
in  a  fashion  that  results  in  groups  of  weights  with  supports  that  are  well 
spaced  in  the  hyperbolic  metric.  One  last  procedure  ensures  that  the  result¬ 
ing  measures  are  symmetric  about  the  real  axis;  this  step  helps  to  minimize 
numerical  errors  that  manifest  themselves  later  as  small  asymmetries  in  the 
final  transfer  functions. 

The  output  from  the  first  program  is  a  file  that  contains  sets  of  ordered 
pairs  of  complex  numbers.  The  sets  correspond  to  the  partition  given  in 
Equations  (26)  and  (27).  The  second  C  program  takes  as  input  the  parti¬ 
tioned  measure  data  from  the  first  Program  and  a  list  of  points  at  which 
to  compute  the  value  of  the  solution.  The  program  constructs  the  Blaschke 
products  and  the  Jones’  interpolating  functions  of  Lemma  6  for  each  part  of 
the  partitioned  measure,  and  evaluates  the  partial  solutions  bk  (z)  at  each  of 
the  data  points.  Summing  the  partial  solutions  at  each  data  point  produces 
the  desired  evaluations  b(z). 

In  summary,  given  a  description  of  the  Carleson  measure,  the  pair  of  C 
programs  will  compute  values  of  approximating  solutions  to  the  inhomoge¬ 
neous  Cauchy  Riemann  equations  on  any  set  of  data  points  in  the  complex 
plane.  The  computed  solutions  converge  in  distribution,  on  any  compact 
set  to  a  bounded  solution  of  the  equations.  The  restriction  of  the  approx- 
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imations  to  the  imaginary  axis  are  uniformly  bounded  L°°  functions  that 
converge  in  the  weak-star  topology  as  the  discretization  of  the  Carleson 
measure  converges  in  distribution  to  the  actual  Carleson  measure. 

4.4  Results 


Figure  5:  Transfer  function  and  pulse  response  for  Q(z)  with  r  =  0.37. 

Figures  5  to  9  show  the  results  of  the  computational  solution  to  the  problem 
described  in  section  2. 

Figure  5  contains  two  graphs  that  describe  the  solution  for  the  H°°  pa¬ 
rameter  Q  in  terms  of  a  transfer  function  and  the  time  domain  response  to  a 
square  input  pulse  of  unit  magnitude  and  1  second  duration.  The  value  for 
the  delay  chosen  was  r  =  0.37.  The  causal  character  of  the  pulse  response 
verifies  that  the  computed  Youla  parameter  is  an  H°°  function.  Figure  6 
contains  Nyquist  plots  of  the  open  loop  transfer  functions  of  the  combined 
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Figure  6:  Nyquist  plot  of  loop  gain  for  controlled  system. 


system  including  Plant  and  feedback  compensator  for  two  values  of  the  delay 
r  =  0.06  and  r  =  0.37.  The  frequency  variable  that  parameterizes  the  curves 
is  y.  A  comparison  with  the  corresponding  Nyquist  plot  from  [EOT92],  Fig¬ 
ure  7,  illustrates  well  the  differences  in  the  controllers  that  are  produced  by 
the  two  different  approaches.  The  controller  of  Enns  et  al.  does  have  better 
stability  and  better  low  frequency  sensitivity,  this  is  to  be  expected,  since 
the  system  that  is  considered  in  the  example  is  sufficiently  simple  that  the 
skew  Toeplitz  theory  from  [OST93]  that  Enns  et  al.  use  is  able  solve  the 
Nehari  problem  with  an  optimal  solution.  In  order  to  do  this  they  require 
that  the  system  be  factored  as  a  product  of  an  H00  function  and  a  rational 
function  with  inverse  in  H00;  further,  each  factor  needs  to  be  decomposed 
by  an  inner  outer  factorization.  Solutions  under  conditions  more  relaxed 
then  those  in  [OST93]  appear  in  the  works  of  Flamm  and  Yang  [FY94]  and 
Ozbay  [093].  Although  a  large  number  of  interesting  systems  satisfy  the 
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conditions  imposed  in  these  works,  the  conditions  are  restrictive.  While 
computation  of  inner-outer  factorizations  is  a  non-trivial  problem,  recent 
work  by  Flamm  and  Crow  [FC94]  which  addresses  the  problem  of  comput¬ 
ing  numerical  approximations  to  inner  outer  factorizations  should  extend 
the  applicability  of  the  results  in  [OST93]  to  more  complicated  examples. 
The  method  presented  in  this  paper  avoids  both  of  these  restrictions  by 
avoiding  the  operator  theoretic  approach  of  [OST93]  entirely. 


Figure  7:  \S(iy)\  =  |(1  +  P{iy)C{iy ))  1 


Figures  7  and  8  plot  the  transfer  functions  that  determine  the  closed 
loop  sensitivity  and  the  robustness.  Comparison  with  Figures  9  and  10  of 
[EOT92] ,  provides  confirmation  of  the  comments  made  in  the  previous  para¬ 
graph.  The  graph  of  the  largest  singular  value  of  P(iy ),  Figure  9,  provides 
an  indication  of  how  far  the  computed  solution  is  from  an  optimal  solution. 
An  optimal  solution  would  have  a  flat  response  with  value  less  than  1. 
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2.5 


x=  0.06 
x-  0.37 


1.5 


Figure  8:  \W2(iy)C(iy)S(iy)\. 


Figure  9:  +  \W2(iy)C(iy)S(iy)\* . 
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5  General  underdetermined  systems  of  linear  Dio- 
phantine  equations 

So  far  this  paper  has  dealt  only  with  controller  design  for  a  single-input 
single-output  plant.  It  turns  out  that  the  same  basic  methods  extend  nat¬ 
urally  to  the  multiple-input  multiple-output  case  at  the  cost  of  an  increase 
in  the  algebraic  complexity  of  the  formulation.  This  section  shows  how  the 
standard,  general  formulation  for  an  H00  control  problem  can  be  reduced 
to  two  general,  underdetermined  systems  of  linear  Diophantine  equations. 
After  extending  the  algebraic  setting  from  Section  3,  Theorem  1  is  used 
to  give  expressions  for  the  solutions  to  these  general  systems  of  equations. 
A  practical  implementation  of  the  methods  of  this  section  would  take  the 
same  form  and  encounter  the  same  problems  as  the  computations  described 
in  Section  4. 


Figure  10:  Configuration  of  Multi-input  Multi-output  Controller 


Figure  10  depicts  a  feedback  controller  for  a  general  class  of  multiple- 
input  multiple-output  robust  control  problems  [Fra87].  The  plant  has  a 
block  transfer  function  matrix 


Fn  F\2 

F21  F22  ’ 


and  the  equations  for  the  closed  loop  system  are: 


z  =  Fnw  +  F12U 
y  —  F21W  +  F22U 
u  =  Gy 

The  robust  stabilization  problem  is  the  problem  of  choosing  a  controller  with 
transfer  function  matrix  G  in  a  way  that  minimizes  the  operator  norm  of 
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the  matrix  transfer  function  that  maps  the  input  signal  w  on  to  the  output 
signal  z.  For  a  stable  linear  system  the  operator  norm  is  the  L00  norm  of  the 
largest  singular  value  of  the  transfer  function  matrix  which,  in  the  single¬ 
input  single-output  case  is  just  the  H°°  norm  of  the  transfer  function.  The 
following  theorem  combines  a  number  of  results  from  Chapter  4  of  [Fra87] 
including  the  Youla  parameterization  of  stabilizing  controllers. 

Theorem  8  [Fra87] 

Assume  that  F  is  stabilizable,  then: 

(i)  G  stabilizes  F  if  and  only  if  G  stabilizes  Fw 

(ii)  Suppose  F22  =  FiF2-1  =  Ff 1  F\  are  co-prime  factorizations  of  F22, 
then  there  exist  X\,  X2,  and  X\,  X2  such  that? 

X2  -Xi  1  [  F2  X\ 

-Pi  F2  \  [  Fi  X2 

and  the  set  of  all  G  stabilizing  F22  is  parameterized  by  the  formulae 

G  =  {-Xl+F2Q){X2-FlQ)~l 

=  (X2-QFlr1(-Xl  +  QF2) 

Q  G  H°°. 

(in)  With  G  given  by  the  parameterization  in  (ii),  and  with  the  transfer 
functions  T\,  T2,  T3  given  by 

T\  =  Fn  —  F12F2X1F21 
T2  —  F12F2 
T3  =  F2F21, 

the  transfer  function  from  w  to  z  in  Figure  10  equals  T\  —  T2QT3. 

Although  this  theorem  was  proved  by  Francis  in  the  setting  of  rational 
H00  functions,  a  suitable  definition  of  co-primeness  broadens  its  scope  to 
more  general  rings.  An  appropriate  notion  of  co-primeness  for  the  ring  H°° 
is  one  based  on  the  condition  in  the  premise  of  Lemma  2;  this  link  between 
the  co-prime  factorizations  for  transfer  functions  and  the  corona  theorem 

51  is  used  to  denote  the  identity  matrix  which,  in  this  context,  is  a  matrix  of  constant 
H°°  functions 
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was  stated  by  Baras  [Bar80].  Also  related  to  the  theorem  is  the  work  of 
Smith  [Smi89]  who  shows  that  any  plant  that  has  a  factorization  over  H00 
and  is  feedback  stabilizable  has  a  doubly  co-prime  factorization. 

Given  a  plant  transfer  function  matrix  F  with  left  and  right  co-prime 
factorizations,  the  design  of  an  optimal  controller  requires  two  steps:  the 
first  is  to  find  eight  H00  matrices  that  satisfy  Equation  (42),  and  the  second 
is  the  search  for  the  H00  matrix  that  minimizes  the  norm  of  the  closed  loop 
transfer  function 

P  =  Ti-  T2QT3.  (43) 

Both  steps  can  be  reformulated  as  the  solution  to  an  underdetermined  sys¬ 
tem  of  equations  of  the  form 

Ax  =  b.  (44) 

Finding  matrices  that  satisfy  Equation  (42)  is  equivalent  to  finding  a  left 
factorization  F  =  F2~lFx  a  right  factorization  F  =  F\  FT1,  and  H°°  matrices 
ATx,  X2,  X\  and  X2  that  satisfy  the  four  equations 


X2F2  —  X\Fx  =  1 

(45) 

r*H 

II 

X 

1 

.  <N 

(46) 

F2Fi  -  F\F2  =  0 

(47) 

o 

II 

1 

* 

. « 

(48) 

Equations  (45)  and  (46)  axe  matrix  Bezout  equations,  and  Equation  (47)  is 
automatically  satisfied  since  the  left  and  right  factorizations  are  factoriza¬ 
tions  of  the  same  transfer  function  matrix.  Given  left  and  right  co-prime 
factorizations  of  F,  and  arbitrary  solutions  Y\  and  Y2  to  Equation  (46) ,  and 
Y\ ,  Y2  to  Equation  (45) ,  a  little  algebraic  manipulation  yields  the  following 
parameterization  of  all  doubly  co-prime  factorizations  that  satisfy  equations 
(45  -  48): 

Xx  =  YX  +  F2A 

X2  =  Y2  +  FxA 

X\  =  Yx  +  (A  —  YxY2  +  Y\Y2)F2 

X2  =  Y2  +  (A  —  YiY2  +  Y{Y2)F2. 

The  parameter  A  is  a  matrix  with  entries  in  H°°.  With  this  result  the 
computation  of  a  doubly  co-prime  factorization  reduces  to  the  solution  of 
the  two  matrix  Bezout  equations  (45)  and  (46).  The  form  of  these  equations 
is  similar  to  the  form  of  (43),  and  with  suitable  substitutions  the  solution 
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of  each  of  the  three  equations  is  subsumed  by  the  following  problem:  given 
A1,  B1,  A2,  B2,  and  C ,  find  X1  and  X2  that  solve 

A1XlBl  +  A2X2B2  =  C.  (49) 

This  equation  is  the  matrix  analog  of  (8)  for  multiple-input  multiple-output 
systems. 

Equation  (49)  has  the  form  of  a  general  linear  equation  in  the  entries 
of  the  matrices  X1  and  X2.  This  equation  is  recast  in  the  form  of  (44), 
Ax  =  b  by  stacking  the  columns  of  the  matrices  X1  and  X2  to  form  a  long 
vector  x ,  stacking  the  columns  of  C  to  form  the  vector  6,  and  replacing 
the  left  and  right  multiplying  matrices  by  one  left  multiplying  matrix  A. 
If  the  dimension  of  A  is  m  x  n  then  m  <  n,  and  A  represents  a  module 
homomorphism  with  domain  H°°x  xH°°,  and  image  H°°x  .*?.  xH°°. 

As  in  Section  3.1  the  solution  to  (44)  is  based  on  Theorem  1,  but  the 
definition  of  the  spaces  and  the  operators  in  (10)  and  Figure  2  need  to  be 
changed.  Define  the  following  modules  over  a  ring  R 

A °{R)  =  Rx  m.  xR 

A.l(R)  =  Rx  ».  xR 

A 2{R)  =  An~m~1(Rx  xR). 

The  three  rings  of  interest  are  the  same  as  those  in  Section  3,  H00,  B,  and 
C.  Denote  the  rows  of  the  matrix  A  by  A\...  Am  and  the  columns  by 
a\  . .  .an,  then  Aj  G  A1(H°°),  and  aj  G  A0(H°°).  Define  the  homomorphism 
Pa  ■  A1(jR)  ->  A°(R)  by  Pax  =  Y^xiai>  then  Equation  (44)  can  be  written 

Pax  =  b.  let  {e^}  be  a  basis  for  An-,n-1  C”,  and  let  y  =  ypep  be  an  element 

of  A 2{R).  Define  a  second  homomorphism  Pa  ■  A 2{R)  — >•  A l(R)  by 

(PAy)  =  ^]^(*(4i  A  •••  A  Am  Aep))  (50) 

P 

in  which  the  star  homomorphism6  is  taken  relative  to  the  normal  Euclidean 
scalar  product  on  Cn.  With  these  definitions  the  sequence  (10)  may  be 
rewritten  as  the  sequence 

A2(R)^4A1(B)^>A°(R)— 40,  (51) 

sThe  star  homomorphism  is  determined  by  its  action  on  homogeneous  forms.  On  these 
it  satisfies  the  equation  *(e;,  A ...  A  e;,. )  A  (e;,  A . . .  A  A . . .  A  e„,  where  the  vectors 

e>  are  the  unit  vectors  in  Cn 
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which  is  exact  at  Al(R).  In  fact  the  sequence  in  (51)  may  be  extended 
leftward  to  form  a  complete  sequence,  but  the  definitions  given  are  already 
enough  for  an  application  of  Theorem  1. 

Theorem  2  provides  an  algorithm  that  produces  a  solution  to  (44)  as 
follows.  First  find  x 1  £  A1  (B)  that  solves  the  equation 

Pax 1  =  b.  (52) 


The  ring  B  is  the  ring  of  distributions  with  boundary  values  in  L°°  that 
was  introduced  in  Section  3.1.  The  solution  needs  to  be  chosen  so  that 
dxl/dz  £  A^C),  where  C  is  the  ring  of  Carleson  measures  that  have  support 
on  the  right  half  plane.  With  this  in  mind  choose  x2  £  A2(C)  to  be  a  solution 
of 

_  2  dz1  /  x 

p*x  =  Si  <M> 

and  x3  £  A 2  (6)  to  be  a  solution  of 


(54) 


It  follows  from  Theorem  1  that  a  solution  to  (44)  is  given  by 


x  =  a;1  —  Pax3- 


(55) 


As  in  the  scalar  case,  the  real  computational  problems  lie  in  inverting 
the  three  operators  d/dz  :  A2  (B)  -4  A 2(C),  Pa  '■  A1  (B)  -4  A0  (B)  and  Pa  : 
A2(C)  -4  AX(C).  Fortunately  though,  the  same  approach  that  was  used  in  the 
scalar  case  may  be  applied  to  systems  of  equations  with  some  extra  algebra. 
The  first  operator,  d/dz  is  the  easiest  to  deal  with,  Equation  (54)  may  be 
solved  by  applying  the  method  of  Section  3.2  to  calculate  each  component 
of  x3  from  the  corresponding  component  of  x2.  The  remaining  operators  are 
inverted  by  using  a  method  due  to  [Rao83]  to  construct  a  left  inverse,  and 
some  algebraic  constructions  that  are  similar  to  those  presented  in  [BS86] 
and  [Str83]. 

First  consider  Equation  (52).  Denote  by  Ay  the  rank  m  minors  of  A , 
then  the  index  7  can  take  one  of  n\/(n—m)\m\  values  that  correspond  to  the 
choices  of  m  columns  from  the  n  columns  of  A.  Provided  that  the  functions 
Ay  satisfy  the  condition  of  Lemma  2,  there  exists  a  partition  of  the  plane 
into  sets  Cly  such  that  if  <j>y  is  the  characteristic  function  of  the  set  fi7,  that 


is, 


4>7{z) 


0  z  G  fi7 

1  z  £  fi7  ’ 
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then  A7  is  bounded  away  from  zero  outside  the  set  Sly.  The  distributional 
derivatives  d(f>y/dz  are  Carleson  measures  supported  on  the  boundaries  <9fly. 
Choose  Gy  =  (f)y/Ay,  then  each  Gy  is  a  bounded  analytic  function  on  the 
interior  of  fly,  is  identically  zero  outside  f 17,  and  has  a  distributional  deriva¬ 
tive  dGy/dz  that  is  a  Carleson  measure  supported  on  the  boundary  dfly. 
Further,  the  functions  Gy  solve  the  equation 


S  Ay  Gy  -  1. 


[Rao83]  uses  the  Cauchy  Binet  theorem  to  show  that  if 


9jk  —  Gy 

7 


dAy 

da,f.j 


then  the  matrix  G  =  [gjk]  is  a  right  inverse  of  A  with  rank  m  minors  Gy.  It 
follows  that  a  solution  to  equation  (52)  is  given  by 


x 


l 

j 


EE 


<^7  9  Ay 
Ay  dakj  k 


(56) 


The  final  equation  that  needs  to  be  solved  is  Equation  (53)  P^x 2  = 
dxl /dz.  Let  x 2  have  components  yp  with  respect  to  the  canonical  basis  for 
A 2(C).  When  the  solution  from  (56)  is  substituted  for  a1,  and  the  expression 
for  the  operator  Pa  from  (50)  is  expanded  in  coordinates,  the  j’th  component 
of  Equation  (53)  becomes 


EvA  = 

P,a 


The  summation  on  the  left  hand  side  in  this  formula  is  taken  over  all  multi¬ 
indices  a  and  such  that  j  £  a  U  /?,  a  fl  (3  =  0,  and  |a|  —  m.  Substituting 
the  solution  for  Gy  gives 


/3, a 


E  E 

k 


bk  dAy  d<py 

Ay  dakj  dz 
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It  follows  from  the  choice  of  (p  that  the  sum  on  the  right  hand  side  is 
supported  entirely  on  the  curve  segments  3fi7p  fl  <9Q7p.  So  the  components 
yp  of  the  solution  x 2  are  measures  supported  on  the  boundaries  <9f27p ,  and  at 
any  point  on  these  boundaries  there  are  n  equations  for  the  n!/m!(n— m— 1)! 
variables  yp  of  the  form 

3  £  Ip  U  7 q 
3  e  7p  -  7 q 

3  e  7p  n  7  q 

The  arbitrary  signs  are  determined  by  the  sense  of  integration  inherent  in 
the  measures  dipp/dz  and  dipq/dz 

Although  the  algebra  associated  with  the  inversion  of  the  operators 
Pa  :  A1  (B)  — >  A0  (13)  and  Pa  ■  A 2(C)  —>  Al(C)  seems  complicated,  the 
real  computational  difficulties  are  the  same  as  those  experienced  with  the 
single-input  single-output  system,  namely,  choosing  a  partition  ft7  and  com¬ 
puting  minimal  norm  solutions  of  db/dz  =  y,  for  a  Carleson  measure  y.  The 
requirement  that  Lemma  2  places  on  the  minors  Ay  of  the  matrix  A  induces 
the  appropriate  co-primeness  conditions  on  the  left  and  right  factorizations 
of  the  transfer  function  matrix  F  for  the  multiple-input  multiple-output 
system. 


yp-A-a  — 


p,0t 


h.  A  A 


1  dK 

Ayq  da,jk 
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6  Conclusion 


This  paper  has  presented  a  new  computational  method  for  H00  controller 
design.  The  method  places  two  requirements  on  the  systems  to  which  it 
applies:  an  explicitly  computable  co-prime  factorization  of  the  system  over 
H00  functions  should  exist,  and  sufficient  information  about  the  location 
of  the  zeros  of  the  factors  is  needed  to  construct  the  partition  of  unity  in 
Lemma  2.  These  requirements  are  very  close  to  necessary  conditions  for 
a  linear  plant  to  be  stabilizable,  a  fact  that  indicates  that  the  techniques 
presented  are  potentially  widely  applicable. 

The  method  has  been  demonstrated  on  a  simple  example  drawn  from  the 
literature.  For  this  example  independent  methods  may  be  used  to  construct 
a  controller  that  is  optimal  in  the  sense  of  H00  control;  for  this  reason  it  pro¬ 
vides  a  good  indication  of  how  close  the  controller  computed  for  the  example 
is  to  an  optimal  controller.  For  a  controller  design  to  be  truly  practical  it  is 
important  that  it  produce  near  optimal  controllers.  This  is  particularly  true 
for  infinite  dimensional  systems  which  often  have  transfer  functions  that  are 
sensitive  to  small  parameter  variations.  Under  this  criterion  for  practicality, 
the  conclusion  to  be  drawn  from  the  result  in  Section  4  is  that  more  work 
is  needed. 

Two  areas  in  need  of  further  work  stand  out.  The  first  concerns  the 
selection  of  the  partition  of  unity  that  is  postulated  in  in  Lemma  2,  and  is 
used  in  the  construction  of  the  Carleson  measures.  The  particular  selection 
made  for  a  given  problem  affects  the  quality  of  the  solution  through  the 
norm  of  the  inverse  in  Equation  (18),  and  through  the  Carleson  constant 
associated  with  the  Blaschke  product  in  the  inequality  (19).  The  intricate 
construction  that  is  required  in  the  Carleson’s  proof  of  the  Corona  Theorem 
would  indicate  that  in  the  most  general  case  choosing  an  optimal  selection 
is  a  difficult  problem,  however,  many  of  the  situations  that  axe  of  interest  in 
engineering  are  described  by  boundary  value  problems  and  delay  differential 
equations  of  the  type  presented  in  this  paper;  in  these  cases  the  additional 
structure  provided  by  the  problem  description  can  be  exploited  to  provide 
partitions  of  unity  without  recourse  to  elaborate  constructions. 

The  second  area  for  future  work  is  the  problem  of  constructing  bounded 
solutions  to  the  Cauchy  Riemann  equations.  In  particular,  attention  should 
be  paid  to  the  interpolating  function  that  is  used  in  the  computation  of 
the  solutions.  Although  the  function  given  in  Theorem  3  is  optimal  in  the 
sense  that  it  satisfies  the  bound  given  in  Equation  (22)  independent  of  the 
measure,  a  particular  choice  of  interpolating  function  tailored  to  a  particular 
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measure  could  produce  a  lower  bound. 

Finally,  a  note  on  the  calculations.  The  method  described  is  computa¬ 
tionally  intensive,  however,  with  careful  programming,  it  certainly  is  fea¬ 
sible.  The  computations  for  the  example  presented  took  minutes,  rather 
than  hours,  on  a  workstation  with  a  RISC  processor  and  floating-point  co¬ 
processor7.  The  predominant  computation  involves  evaluating  a  small  set 
of  functions  over  a  large  set  of  data  points  with  no  interdependencies  in  the 
evaluations.  This  type  of  calculation  is  trivially  parallelizable  on  massively 
parallel  architectures. 


7The  calculations  were  made  in  compiled  C  code  on  a  SPARC  Station  10 
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